function y=counterfactual_stderr(beta, Scenario,mc, counterfactual)

N=size(beta,2);
S=size(beta,3);
L=size(beta,1);
option=L-size(Scenario,2)+1;

% z=zeros(option,2,S-5000);
z=zeros(option,2,size(1:50:S,2));
% Scenario(1,1)=fsolve('price',4,...
%     [optimset('TolFun',1e-5, 'MaxFunEvals', 100000)],...
%     mc,beta,Scenario, counterfactual);
%Scenario(1,1)
w=[eye(option-1) Scenario];


w(option,:)=0;

G=1;
for s=1:G:S
    share=zeros(option,1);
    
    for n=1:N
        numer=exp(w*beta(:,n,s));
        numer=numer.*[counterfactual 1]';
        share=share+numer/sum(numer);
    end
    
% z(:,1,s)=share/(N);
% 
% z(:,2,s)=(Scenario(1,1)-mc)*z(:,1,s);

z(:,1,1+(s-1)/G)=share/(N);

z(:,2,1+(s-1)/G)=(Scenario(1,1)-mc)*z(:,1,1+(s-1)/G);
end

y=reshape(z,option*2,[])';

